Spatial Analysis in R

2020-02-19

Learning Objectives

Plot vector geospatial data (points/polygons)
Plot raster geospatial data (image)
Extracting data from raster
Making maps (Compass rose, overlays, etc)

Motivation

We are researchers who want to answer the question "What range of environmental conditions do taxon tend to inhabit?"

General workflow:

Plot locations of taxon
Plot map of environmental conditions
Plot contextual information
Extract data about environmental conditions at those locations
Variants of this question/workflow are applicable to others study systems.

–> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –>

Representation of spatial data

Vector data model uses points, lines, and polygons (e.g. locations of animals, ecosystem types)
Raster data model divides surface into cells of constant size (e.g. gridded temperature)
Non spatial attribute data
How many counts at a point?
Length of a transect
Land-cover type of a polygon

To answer our research question, we will need to:

Read/write spatial data
Represent geographic and attribute data
Transform between different models of Earth
Geographic operations
Make maps

Spatial analysis in R

R packages for spatial analysis

Many tools for spatial analysis in R
rgdal released in 2003–import from more geographic data formats
sp released in 2005–creates spatial objects with classes and generic methods for points, lines, polygons, grids, and attribute data (but still lacked ability to do geometric operations)
raster released in 2010–support for raster operations and more

R packages linking to GIS software

GRASS, spgrass6, and rgrass7
RSAGA, RPyGeo, and RQGIS

Visualization

ggplot2
ggmap
rasterVis
tmap
leaflet
mapview

My choice: `sf`

Developed in 2016 with support from R consortium
Based on simple features ISO standards that are not R specific
Replaces sp, rdgal for read/write, and rgeos for geometric operations
Integrates well with tidyverse

Resources

Simple Features package
Vignettes, blogs, presentations
Cheatsheet
issue tracker
r-spatial Edzer Pebesma's blog
Spatial Data Science
Geocomputation with R
excellent in-depth treatment
includes theory, examples, and covers rasters

PART I: PLOTTING POINT AND SHAPEFILE DATA

DATA: 1. taxon point locations, 2. relevant shapefile (Ecosystem type?, country?)

`sf` package

Simple features

feature "abstraction of real world phenomena"
simple feature "feature with all geometric attributes described by piecewise straight lines or planar interpolation between sets of points"
Represents geometry using points, lines, polygons, or collections of those features
ISO standard supported by lots of software

Basic structure

sf (simple feature) objects are extended data.frames or tibbles
attribute data stored as tibble
geometry stored as a list column
can have multiple types of geometry in one object
class is sfc (simple feature columns) with bounding box bbox and CRS (coordinate reference system) attributes

## Reading layer `NAM-level_1' from data source `/Users/michelebuonanduci/Desktop/QERM597-making-maps/data/NAM-level_1.shp' using driver `ESRI Shapefile'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.

## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Value 'mp\zip\NAM-1' of field NAM-level_1.PERIMETER parsed
## incompletely to real 0.
## Simple feature collection with 19 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 11.71563 ymin: -28.96946 xmax: 25.2567 ymax: -16.95989
## epsg (SRID):    NA
## proj4string:    NA
## [1] "sf"         "data.frame"
## Simple feature collection with 6 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 13.62839 ymin: -28.96946 xmax: 25.2567 ymax: -17.47639
## epsg (SRID):    NA
## proj4string:    NA
##   AREA PERIMETER      ID CAPTION                       geometry
## 1   NA         0 Caprivi Caprivi POLYGON ((24.77467 -17.8854...
## 2   NA         0  Erongo  Erongo POLYGON ((15.81059 -23.2945...
## 3   NA         0  Hardap  Hardap POLYGON ((18.48226 -25.4376...
## 4   NA         0   Karas   Karas POLYGON ((16.3988 -25.64753...
## 5   NA         0   Karas    <NA> POLYGON ((15.20197 -27.0306...
## 6   NA         0   Karas    <NA> POLYGON ((14.92512 -26.3011...
## Reading layer `enproads' from data source `/Users/michelebuonanduci/Desktop/QERM597-making-maps/data/enproads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 565 features and 3 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008
## epsg (SRID):    NA
## proj4string:    NA
## [1] "sf"         "data.frame"
## Simple feature collection with 6 features and 3 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 14.46841 ymin: -19.41752 xmax: 14.58964 ymax: -19.33077
## epsg (SRID):    NA
## proj4string:    NA
##     TYPE LENGTH KM                       geometry
## 1  Track    0.7  1 LINESTRING (14.50502 -19.35...
## 2  Track    3.6  4 LINESTRING (14.49077 -19.36...
## 3 Graded    2.3  2 LINESTRING (14.4864 -19.417...
## 4  Track    7.1  7 LINESTRING (14.58498 -19.34...
## 5  Track    1.2  1 LINESTRING (14.58498 -19.34...
## 6 Graded    0.1  0 LINESTRING (14.52596 -19.33...

`sf` and `tidyverse`

sf functions begin with st_
Methods for summary, plot
sf methods for filter, arrange, distinct, group_by, ungroup, mutate

methods(class = "sf")
##   [1] [                          [[<-                      
##   [3] $<-                        aggregate                 
##   [5] anti_join                  arrange                   
##   [7] as.data.frame              cbind                     
##   [9] coerce                     dbDataType                
##  [11] dbWriteTable               df_spatial                
##  [13] distinct                   extent                    
##  [15] extract                    filter                    
##  [17] full_join                  gather                    
##  [19] group_by                   group_map                 
##  [21] group_split                identify                  
##  [23] initialize                 inner_join                
##  [25] left_join                  mask                      
##  [27] merge                      mutate                    
##  [29] nest                       plot                      
##  [31] print                      raster                    
##  [33] rasterize                  rbind                     
##  [35] rename                     right_join                
##  [37] sample_frac                sample_n                  
##  [39] select                     semi_join                 
##  [41] separate                   show                      
##  [43] slice                      slotsFromS3               
##  [45] spread                     st_agr                    
##  [47] st_agr<-                   st_area                   
##  [49] st_as_sf                   st_bbox                   
##  [51] st_boundary                st_buffer                 
##  [53] st_cast                    st_centroid               
##  [55] st_collection_extract      st_convex_hull            
##  [57] st_coordinates             st_crop                   
##  [59] st_crs                     st_crs<-                  
##  [61] st_difference              st_force_polygon_cw       
##  [63] st_geometry                st_geometry<-             
##  [65] st_interpolate_aw          st_intersection           
##  [67] st_intersects              st_is_polygon_cw          
##  [69] st_is                      st_line_merge             
##  [71] st_linesubstring           st_make_valid             
##  [73] st_minimum_bounding_circle st_nearest_points         
##  [75] st_node                    st_normalize              
##  [77] st_point_on_surface        st_polygonize             
##  [79] st_precision               st_segmentize             
##  [81] st_set_precision           st_simplify               
##  [83] st_snap_to_grid            st_snap                   
##  [85] st_split                   st_subdivide              
##  [87] st_sym_difference          st_transform_proj         
##  [89] st_transform               st_triangulate            
##  [91] st_union                   st_voronoi                
##  [93] st_wrap_dateline           st_write                  
##  [95] st_zm                      summarise                 
##  [97] transmute                  ungroup                   
##  [99] unite                      unnest                    
## see '?methods' for accessing help and source code

st_geometry(Namibia) %>% plot

plot(st_geometry(roads))

# Default methods for objects
plot(Namibia)
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

plot(Namibia[,3])

plot(roads)

plot(roads[,1])

plot(filter(roads[,3], KM > 20))

Reading and writing data

Reading shapefiles

fileloc <- system.file("shape/nc.shp", package = "sf")
nc <- read_sf(fileloc)

st_crs() used to convert coordinates (do projections) and do datum transformations
Proj.4 string or ESPG code contains this info
Web repositories for ESPG and Proj.4 information

st_crs("+proj=longlat +datum=WGS84") # Proj.4 string"
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
st_crs(3857)                         # ESPG code: public registry of spatial reference systems, 3857 = web-based mapping (i.e. Google                                         maps, OpenStreetMap)
## Coordinate Reference System:
##   EPSG: 3857 
##   proj4string: "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
st_crs(3857)$units                   # check units
## [1] "m"
st_crs(NA)                           # unknown (assumed planar/Cartesian)
## Coordinate Reference System: NA

units::set_units(a1, km^2)
## 1137.389 [km^2]
units::set_units(a2, km^2)
## 1137.577 [km^2]
units::set_units(a3, km^2)
## 1137.577 [km^2]

Other data sources

Convert from other Spatial* objects using st_as_sf
Convert long/lat data
Packages such as rnaturalearth

Convert long/lat data

head(hotsprings_df)
## # A tibble: 6 x 8
##   STATE   LAT  LONG SpringName          Temp_F Temp_C AREA    USGS_quad    
##   <chr> <dbl> <dbl> <chr>                <dbl>  <dbl> <chr>   <chr>        
## 1 CA     38.8 -123. LITTLE        GEYS…    210     99 SANTA … (WHISPERING …
## 2 CA     41.5 -120. HOT        SPRINGS…    208     98 ALTURAS CEDARVILLE 15
## 3 CA     36.0 -118. COSO        HOT SP…    207     97 DEATH … HAIWEE RESER…
## 4 CA     41.7 -120. LAKE        CITY H…    207     97 ALTURAS CEDARVILLE 15
## 5 CA     36.0 -118. DEVILS        KITC…    207     97 DEATH … HAIWEE RESER…
## 6 CA     40.4 -121. TERMINAL        GE…    205     96 SUSANV… MT. HARKNESS…

Convert long/lat data

Assume that these data are WGS84 (crs = 4326)
WGS84: World Geodic System 1984
4326=coordinate system based on Earth's center of mass

hotsprings <- st_as_sf(hotsprings_df,
                       coords = c("LONG", "LAT"),
                       crs = 4326)
head(hotsprings)
## Simple feature collection with 6 features and 6 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -122.748 ymin: 36.036 xmax: -117.769 ymax: 41.67
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 6 x 7
##   STATE SpringName  Temp_F Temp_C AREA   USGS_quad                 geometry
##   <chr> <chr>        <dbl>  <dbl> <chr>  <chr>                  <POINT [°]>
## 1 CA    LITTLE    …    210     99 SANTA… (WHISPER…        (-122.748 38.767)
## 2 CA    HOT       …    208     98 ALTUR… CEDARVIL…        (-120.078 41.534)
## 3 CA    COSO      …    207     97 DEATH… HAIWEE R…        (-117.769 36.045)
## 4 CA    LAKE      …    207     97 ALTUR… CEDARVIL…         (-120.206 41.67)
## 5 CA    DEVILS    …    207     97 DEATH… HAIWEE R…        (-117.802 36.036)
## 6 CA    TERMINAL  …    205     96 SUSAN… MT. HARK…        (-121.375 40.421)

Data from packages

world <- ne_countries(scale = "medium", returnclass = "sf")
st_geometry(world) %>% plot()

Writing data

shapefiles
database connections
use st_write

Exercises

Load the roads shapefile (enproads.shp) and upload Zebra.csv and plot the geometry.
If we look at our roads data, we see that there are several types of roads in the shapefile. How close or far from different road types do zebra move?

Exercise 1

roads <- st_read(here("data", "enproads.shp"), crs = "+init=epsg:4326") %>% #4326= coordinate system based on Earth's center of mass
  st_transform("+init=epsg:32733") #32733 = spatial reference for Nambia 
## Reading layer `enproads' from data source `/Users/michelebuonanduci/Desktop/QERM597-making-maps/data/enproads.shp' using driver `ESRI Shapefile'
## Simple feature collection with 565 features and 3 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 14.35337 ymin: -19.48418 xmax: 17.15124 ymax: -18.5008
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

st_geometry(roads) %>% plot


zebra_sf <- read.csv(here("data", "Zebra.csv")) %>% 
  dplyr::select(ID = Name, 4:6) %>% 
  mutate(timestamp = as.POSIXct(lubridate::mdy_hm(Date))) %>%
  st_as_sf(., coords = 3:4, crs = "+init=epsg:4326") %>% #longlat, converting foreign object to SF object
  st_transform("+init=epsg:32733") #convert to UTM

ggplot() +
  geom_sf(data=roads) +
  geom_sf(data=zebra_sf, aes(color=ID))

Exercise 2

How close or far from different road types do zebra move?

unique(roads$TYPE)
## [1] Track  Graded Gravel Tar   
## Levels: Graded Gravel Tar Track

#Filter roads based off type:
large_roads <- filter(roads, TYPE %in% c("Tar", "Gravel"))
small_roads <- filter(roads, TYPE %in% c("Graded", "Track"))

ggplot() +
  geom_sf(data=large_roads, size=1.5) + 
  geom_sf(data=small_roads, size=0.6) + 
  geom_sf(data=zebra_sf, aes(color=ID))


# Find the minimum distance (in meters) of each point to a large road
large_dist<- st_distance(y=zebra_sf, x=large_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds
zebra_sf$large_road_dist <- apply(large_dist, 2, min)

# Find the minimum distance (in meters) of each point to a small road
small_dist<- st_distance(y=zebra_sf, x=small_roads) # a units matrix dim =[nrow(x), nrow(y)]; takes 10-20 seconds
zebra_sf$small_road_dist <- apply(small_dist, 2, min)

head(data.frame(zebra_sf))
##      ID            Date           timestamp                 geometry
## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501)
## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00   POINT (596576 7879266)
## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260)
## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00   POINT (584733 7890124)
## 5 AG059  4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269)
## 6 AG059  4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266)
##   large_road_dist small_road_dist
## 1        7.845860       3479.8987
## 2      163.187201        241.8660
## 3      156.746071        245.8446
## 4        9.115608       1605.2154
## 5      151.769751        227.9759
## 6      151.358102        231.6984

ggplot(zebra_sf) +
  geom_histogram(aes(large_road_dist, fill="large roads"), alpha=0.5) +
  geom_histogram(aes(small_road_dist, fill="small roads"), alpha=0.5) +
  scale_fill_manual(labels=c("large roads", "small roads"), values=c("blue", "orange"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

PART II: PLOTTING RASTER DATA AND EXTRACTING INFORMATION

DATA: 1. taxon point locations, 2. relevant shapefile (Ecosystem type?, country?)

Reading raster data

## class      : RasterLayer 
## dimensions : 474, 1274, 603876  (nrow, ncol, ncell)
## resolution : 231.6564, 231.6564  (x, y)
## extent     : 431873.8, 727004, 7844489, 7954294  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /Users/michelebuonanduci/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif 
## names      : ndvi_mean_utm 
## values     : -558.117, 5303.021  (min, max)
## class      : RasterLayer 
## dimensions : 474, 1274, 603876  (nrow, ncol, ncell)
## resolution : 231.6564, 231.6564  (x, y)
## extent     : 431873.8, 727004, 7844489, 7954294  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /Users/michelebuonanduci/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif 
## names      : ndvi_mean_utm 
## values     : -558.117, 5303.021  (min, max)
## class      : RasterLayer 
## dimensions : 474, 1274, 603876  (nrow, ncol, ncell)
## resolution : 231.6564, 231.6564  (x, y)
## extent     : 431873.8, 727004, 7844489, 7954294  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /Users/michelebuonanduci/Desktop/QERM597-making-maps/data/ndvi_mean_utm.tif 
## names      : ndvi_mean_utm 
## values     : -558.117, 5303.021  (min, max)

Converting projection

Writing raster data

Extracting data

EXERCISES

Load the NDVI raster file and plot. Overlay Zebra location
Extract NDVI at zebra locations
Plot zebra location data. Color points according to NDVI values

Exercise 1

## class      : RasterLayer 
## dimensions : 474, 1274, 603876  (nrow, ncol, ncell)
## resolution : 231.6564, 231.6564  (x, y)
## extent     : 431873.8, 727004, 7844489, 7954294  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=33 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : ndvi_mean_utm 
## values     : -0.0558117, 0.5303021  (min, max)

Exercise 2

##      ID            Date           timestamp                 geometry
## 1 AG059 4/20/2009 15:40 2009-04-20 15:40:00 POINT (594243.6 7885501)
## 2 AG059 4/20/2009 13:01 2009-04-20 13:01:00   POINT (596576 7879266)
## 3 AG059 4/20/2009 12:02 2009-04-20 12:02:00 POINT (596579.5 7879260)
## 4 AG059 4/20/2009 17:00 2009-04-20 17:00:00   POINT (584733 7890124)
## 5 AG059  4/21/2009 1:40 2009-04-21 01:40:00 POINT (596611.3 7879269)
## 6 AG059  4/21/2009 1:00 2009-04-21 01:00:00 POINT (596605.8 7879266)
##   large_road_dist small_road_dist NDVI_mean
## 1        7.845860       3479.8987 0.2409661
## 2      163.187201        241.8660 0.2682169
## 3      156.746071        245.8446 0.2682169
## 4        9.115608       1605.2154 0.2380762
## 5      151.769751        227.9759 0.2419201
## 6      151.358102        231.6984 0.2419201

Acknowledge existance of multidimenional data and packages. netCDF files

PART III: ADDING BASE MAPS

Adding basemaps!

Focus on ggmap
As of mid-2018, requires registering with Google and obtaining an API key
Requires providing a valid credit card (yikes!), though charges are nonexistent or at least very minimal
See the ggmap GitHub page for more information about API keys

Basemap options from `ggmap`

Basemaps available from Google, Stamen, or Open Street Map
Terrain, satellite, or watercolor
See this helpful cheat sheet

Making maps using `ggmap`

Two steps:

Download basemap raster
Plot raster and overlay other spatial data

Register API key at start of session

This is my personal API key, which you can use for the purposes of this class!

register_google(key = "AIzaSyBRkw_wzDKuWZDikdD86Kmp1Sa6PzuCKFc")

Stamen basemaps

To add a Stamen basemap, first define the bounding box for the basemap you would like to download.

zebra_bbox <- c(14, -20, 17.5, -18)
names(zebra_bbox) <- c("left", "bottom", "right", "top")

Stamen basemaps: terrain

terr_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="terrain")

ggmap(terr_basemap) +
  geom_sf(data=roads, inherit.aes = FALSE) + 
  geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
  coord_sf(crs = st_crs(4326))

Stamen basemaps: watercolor

wat_basemap <- get_stamenmap(bbox=zebra_bbox, zoom=9, maptype="watercolor")

ggmap(wat_basemap) +
  geom_sf(data=roads, inherit.aes = FALSE) + 
  geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
  coord_sf(crs = st_crs(4326))

Google basemaps

To add a Google basemap, first define the center coordinates for the basemap you would like to download.

zebra_center <- c(15.8, -19)
names(zebra_center) <- c("lon", "lat")

Google basemaps: satellite

sat_basemap <- get_googlemap(center=zebra_center, zoom=8, size=c(640, 640), scale=2,
                             maptype="satellite")
ggmap(sat_basemap) +
  geom_sf(data=roads, inherit.aes = FALSE) + 
  geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
  coord_sf(crs = st_crs(4326))

Google basemaps: satellite

Notice that get_googlemap() returns a square basemap tile, which then sets the plotting limits of your map.

To manually set different plotting limits, pass in an empty "base layer" and set x and y limits using ggplot().

Google basemaps: satellite

ggmap(sat_basemap, maprange=FALSE, base_layer=ggplot(aes(x=1, y=1), data=NULL)) +
  xlim(14.3, 17.4) + ylim(-20, -18) + xlab("lon") + ylab("lat") +
  geom_sf(data=roads, inherit.aes = FALSE) + 
  geom_sf(data=zebra_sf, aes(color=ID), inherit.aes = FALSE) +
  coord_sf(crs = st_crs(4326))

Adding basemaps: alternative approaches

An alternative package is RgoogleMaps
No need to get Google API key (yay!)
Does not seems to be compatible with ggplot2 (boo)
Another approach is to download raster data directly from Natural Earth
Have to download manually from website and import tif file as a raster brick, then display red/green/blue values
Only works well if you are mapping a large spatial extent (lacks fine resolution)

PART IV: Making it pretty for publication

Things to tweak

Fix the legend
Add a scale bar and north arrow
Annotations
Inset maps

Natural earth package to access other data

## OGR data source with driver: ESRI Shapefile 
## Source: "/private/var/folders/vn/q1jnls511r35nvbxkfgqzn4h0000gn/T/RtmpqXdESv", layer: "ne_10m_rivers_lake_centerlines"
## with 1455 features
## It has 34 fields
## Integer64 fields read as strings:  rivernum ne_id
## Warning in rgdal::readOGR(destdir, file_name, encoding = "UTF-8",
## stringsAsFactors = FALSE, : Dropping null geometries: 1453
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant
## throughout all geometries

Learning Objectives

Motivation

Representation of spatial data

To answer our research question, we will need to:

Spatial analysis in R

R packages for spatial analysis

R packages linking to GIS software

Visualization

My choice: sf

Resources

PART I: PLOTTING POINT AND SHAPEFILE DATA

DATA: 1. taxon point locations, 2. relevant shapefile (Ecosystem type?, country?)

sf package

Simple features

Basic structure

sf and tidyverse

Reading and writing data

Reading shapefiles

Other data sources

Convert long/lat data

Convert long/lat data

Data from packages

Writing data

Exercises

Exercise 1

Exercise 2

PART II: PLOTTING RASTER DATA AND EXTRACTING INFORMATION

DATA: 1. taxon point locations, 2. relevant shapefile (Ecosystem type?, country?)

Reading raster data

Converting projection

Writing raster data

Extracting data

EXERCISES

Exercise 1

Exercise 2

Acknowledge existance of multidimenional data and packages. netCDF files

PART III: ADDING BASE MAPS

Adding basemaps!

Basemap options from ggmap

Making maps using ggmap

Register API key at start of session

Stamen basemaps

Stamen basemaps: terrain

Stamen basemaps: watercolor

Google basemaps

Google basemaps: satellite

Google basemaps: satellite

Google basemaps: satellite

Adding basemaps: alternative approaches

PART IV: Making it pretty for publication

Things to tweak

Natural earth package to access other data

Fix the legend

Add scale bar and north arrow

Add annotation

My choice: `sf`

`sf` package

`sf` and `tidyverse`

Basemap options from `ggmap`

Making maps using `ggmap`